      subroutine INT1_coord_transform(xi,et,ze,xl,xs,xsi)
!
!     shape functions and derivatives for a 20-node quadratic
!     isoparametric brick element. Watch out: Was 0<=xi<=2, -1<=et,ze<=1 !!!
!                                             Now -1<=xi,et,ze<=1
!
      implicit none
!
      integer i,j,k
!
      real*8 shp(4,20),xs(3,3),xsi(3,3),xl(3,20)
!
      real*8 xi,et,ze,xsj

      real*8 row_norm
!
!     shape functions and their glocal derivatives
!
!     local derivatives of the shape functions: xi-derivative
!
      xi = xi + 1               ! shift the centered shape function
      shp(1, 1)=-0.125+et*(0.250+et*(-0.125-ze*0.125)
     &  +ze*(0.375+ze*0.125))+ze*(-0.250-ze*0.125)
     &  +xi*(0.250+et*(-0.250-0.250*ze)+ze*0.250)
      shp(1, 2)=-0.375+et*(0.250+et*(0.125+ze*0.125)
     &  +ze*(0.125-ze*0.125))+ze*(-0.250+ze*0.125)
     &  +xi*(0.250+et*(-0.250-0.250*ze)+ze*0.250)
      shp(1, 3)=-0.375+et*(0.250+et*(0.125-ze*0.125)
     &  +ze*(-0.125-ze*0.125))+ze*(0.250+ze*0.125)
     &  +xi*(0.250+et*(-0.250+ze*0.250)-0.250*ze)
      shp(1, 4)=-0.125+et*(0.250+et*(-0.125+ze*0.125)
     &  +ze*(-0.375+ze*0.125))+ze*(0.250-ze*0.125)
     &  +xi*(0.250+et*(-0.250+ze*0.250)-0.250*ze)
      shp(1, 5)=-0.125+et*(-0.250+et*(-0.125-ze*0.125)
     &  +ze*(-0.375-ze*0.125))+ze*(-0.250-ze*0.125)
     &  +xi*(0.250+et*(0.250+ze*0.250)+ze*0.250)
      shp(1, 6)=-0.375+et*(-0.250+et*(0.125+ze*0.125)
     &  +ze*(-0.125+ze*0.125))+ze*(-0.250+ze*0.125)
     &  +xi*(0.250+et*(0.250+ze*0.250)+ze*0.250)
      shp(1, 7)=-0.375+et*(-0.250+et*(0.125-ze*0.125)
     &  +ze*(0.125+ze*0.125))+ze*(0.250+ze*0.125)
     &  +xi*(0.250+et*(0.250-0.250*ze)-0.250*ze)
      shp(1, 8)=-0.125+et*(-0.250+et*(-0.125+ze*0.125)
     &  +ze*(0.375-ze*0.125))+ze*(0.250-ze*0.125)
     &  +xi*(0.250+et*(0.250-0.250*ze)-0.250*ze)
      shp(1, 9)= 0.500+et*(-0.500-0.500*ze)+ze*0.500
     &  +xi*(-0.500+et*(0.500+ze*0.500)-0.500*ze)
      shp(1,10)= 0.250+et*(-0.250+ze*(ze*0.250))+ze*(-ze*0.250)
      shp(1,11)= 0.500+et*(-0.500+ze*0.500)-0.500*ze
     &  +xi*(-0.500+et*(0.500-0.500*ze)+ze*0.500)
      shp(1,12)=-0.250+et*(0.250+ze*(-ze*0.250))+ze*(ze*0.250)
      shp(1,13)= 0.500+et*(0.500+ze*0.500)+ze*0.500
     &  +xi*(-0.500+et*(-0.500-0.500*ze)-0.500*ze)
      shp(1,14)= 0.250+et*(0.250+ze*(-ze*0.250))+ze*(-ze*0.250)
      shp(1,15)= 0.500+et*(0.500-0.500*ze)-0.500*ze
     &  +xi*(-0.500+et*(-0.500+ze*0.500)+ze*0.500)
      shp(1,16)=-0.250+et*(-0.250+ze*(ze*0.250))+ze*(ze*0.250)
      shp(1,17)=-0.250+et*(et*(0.250+ze*0.250))-0.250*ze
      shp(1,18)= 0.250+et*(et*(-0.250-ze*0.250))+ze*0.250
      shp(1,19)= 0.250+et*(et*(-0.250+ze*0.250))-0.250*ze
      shp(1,20)=-0.250+et*(et*(0.250-ze*0.250))+ze*0.250
!
!     local derivatives of the shape functions: eta-derivative
!
      shp(2, 1)=et*(0.500
     &  +ze*0.500)+ze*(-0.250-ze*0.250)+xi*(0.250+et*(-0.250
     &  -0.250*ze)+ze*(0.375+ze*0.125)+xi*(-0.125-0.125*ze))
      shp(2, 2)=xi*(0.250+et*(0.250+ze*0.250)+ze*(0.125-ze*0.125)
     &  +xi*(-0.125-0.125*ze))
      shp(2, 3)=xi*(0.250+et*(0.250-0.250*ze)+ze*(-0.125-ze*0.125)
     &  +xi*(-0.125+ze*0.125))
      shp(2, 4)=et*(0.500-0.500*ze)+ze*(0.250-ze*0.250)
     &  +xi*(0.250+et*(-0.250+ze*0.250)+ze*(-0.375+ze*0.125)
     &  +xi*(-0.125+ze*0.125))
      shp(2, 5)=et*(0.500+ze*0.500)+ze*(0.250+ze*0.250)
     &  +xi*(-0.250+et*(-0.250-0.250*ze)+ze*(-0.375-ze*0.125)
     &  +xi*(0.125+ze*0.125))
      shp(2, 6)=xi*(-0.250+et*(0.250+ze*0.250)+ze*(-0.125+ze*0.125)
     &  +xi*(0.125+ze*0.125))
      shp(2, 7)=xi*(-0.250+et*(0.250-0.250*ze)+ze*(0.125+ze*0.125)
     &  +xi*(0.125-0.125*ze))
      shp(2, 8)=et*(0.500-0.500*ze)+ze*(-0.250+ze*0.250)
     &  +xi*(-0.250+et*(-0.250+ze*0.250)+ze*(0.375-ze*0.125)
     &  +xi*(0.125-0.125*ze))
      shp(2, 9)=xi*(-0.500-0.500*ze+xi*(0.250+ze*0.250))
      shp(2,10)=xi*(-0.250+ze*(ze*0.250))
      shp(2,11)=xi*(-0.500+ze*0.500+xi*(0.250-0.250*ze))
      shp(2,12)=-0.500+ze*(ze*0.500)+xi*(0.250+ze*(-ze*0.250))
      shp(2,13)=xi*(0.500+ze*0.500+xi*(-0.250-0.250*ze))
      shp(2,14)=xi*(0.250+ze*(-ze*0.250))
      shp(2,15)=xi*(0.500-0.500*ze+xi*(-0.250+ze*0.250))
      shp(2,16)= 0.500+ze*(-ze*0.500)+xi*(-0.250+ze*(ze*0.250))
      shp(2,17)=et*(-1.000+ze*(-1.000))+xi*(et*(0.500+ze*0.500))
      shp(2,18)=xi*(et*(-0.500-0.500*ze))
      shp(2,19)=xi*(et*(-0.500+ze*0.500))
      shp(2,20)=et*(-1.000+ze*( 1.000))+xi*(et*(0.500-0.500*ze))
!
!     local derivatives of the shape functions: zeta-derivative
!
      shp(3, 1)=et*(-0.250+et*0.250-0.500*ze)+ze*0.500
     &  +xi*(-0.250+et*(0.375-0.125*et+ze*0.250)-0.250*ze
     &  +xi*(0.125+et*(-0.125)))
      shp(3, 2)=xi*(-0.250+et*(0.125+et*0.125-0.250*ze)+ze*0.250
     &  +xi*(0.125+et*(-0.125)))
      shp(3, 3)=xi*(0.250+et*(-0.125-0.125*et-0.250*ze)+ze*0.250
     &  +xi*(-0.125+et*(0.125)))
      shp(3, 4)=et*(0.250-0.250*et-0.500*ze)+ze*0.500
     &  +xi*(0.250+et*(-0.375+et*0.125+ze*0.250)-0.250*ze
     &  +xi*(-0.125+et*(0.125)))
      shp(3, 5)=et*(0.250+et*0.250+ze*0.500)+ze*0.500
     &  +xi*(-0.250+et*(-0.375-0.125*et-0.250*ze)-0.250*ze
     &  +xi*(0.125+et*(0.125)))
      shp(3, 6)=xi*(-0.250+et*(-0.125+et*0.125+ze*0.250)+ze*0.250
     &  +xi*(0.125+et*(0.125)))
      shp(3, 7)=xi*(0.250+et*(0.125-0.125*et+ze*0.250)+ze*0.250
     &  +xi*(-0.125+et*(-0.125)))
      shp(3, 8)=et*(-0.250-0.250*et+ze*0.500)+ze*0.500
     &  +xi*(0.250+et*(0.375+et*0.125-0.250*ze)-0.250*ze
     &  +xi*(-0.125+et*(-0.125)))
      shp(3, 9)=xi*(0.500+et*(-0.500)+xi*(-0.250+et*(0.250)))
      shp(3,10)=xi*(et*(ze*0.500)-0.500*ze)
      shp(3,11)=xi*(-0.500+et*(0.500)+xi*(0.250+et*(-0.250)))
      shp(3,12)=et*(ze*( 1.000))+ze*(-1.000)+xi*(et*(-0.500*ze)+
     &  ze*0.500)
      shp(3,13)=xi*(0.500+et*(0.500)+xi*(-0.250+et*(-0.250)))
      shp(3,14)=xi*(et*(-0.500*ze)-0.500*ze)
      shp(3,15)=xi*(-0.500+et*(-0.500)+xi*(0.250+et*(0.250)))
      shp(3,16)=et*(ze*(-1.000))+ze*(-1.000)+xi*(et*(ze*0.500)+ze*0.500)
      shp(3,17)= 0.500+et*(et*(-0.500))+xi*(-0.250+et*(et*0.250))
      shp(3,18)=xi*(0.250+et*(et*(-0.250)))
      shp(3,19)=xi*(-0.250+et*(et*0.250))
      shp(3,20)=-0.500+et*(et*0.500)+xi*(0.250+et*(et*(-0.250)))
c$$$!
c$$$!     shape functions
c$$$!
c$$$      shp(4, 1)=-0.250+et*(et*(0.250+ze*0.250)
c$$$     &  +ze*(-0.250-ze*0.250))+ze*(ze*0.250)
c$$$     &  +xi*(-0.125+et*(0.250+et*(-0.125-ze*0.125)
c$$$     &  +ze*(0.375+ze*0.125))+ze*(-0.250-ze*0.125)
c$$$     &  +xi*(0.125+et*(-0.125-0.125*ze)+ze*0.125))
c$$$      shp(4, 2)=xi*(-0.375+et*(0.250+et*(0.125+ze*0.125)
c$$$     &  +ze*(0.125-ze*0.125))+ze*(-0.250+ze*0.125)
c$$$     &  +xi*(0.125+et*(-0.125-0.125*ze)+ze*0.125))
c$$$      shp(4, 3)=xi*(-0.375+et*(0.250+et*(0.125-ze*0.125)
c$$$     &  +ze*(-0.125-ze*0.125))+ze*(0.250+ze*0.125)
c$$$     &  +xi*(0.125+et*(-0.125+ze*0.125)-0.125*ze))
c$$$      shp(4, 4)=-0.250+et*(et*(0.250-ze*0.250)
c$$$     &  +ze*(0.250-ze*0.250))+ze*(ze*0.250)
c$$$     &  +xi*(-0.125+et*(0.250+et*(-0.125+ze*0.125)
c$$$     &  +ze*(-0.375+ze*0.125))+ze*(0.250-ze*0.125)
c$$$     &  +xi*(0.125+et*(-0.125+ze*0.125)-0.125*ze))
c$$$      shp(4, 5)=-0.250+et*(et*(0.250+ze*0.250)
c$$$     &  +ze*(0.250+ze*0.250))+ze*(ze*0.250)
c$$$     &  +xi*(-0.125+et*(-0.250+et*(-0.125-ze*0.125)
c$$$     &  +ze*(-0.375-ze*0.125))+ze*(-0.250-ze*0.125)
c$$$     &  +xi*(0.125+et*(0.125+ze*0.125)+ze*0.125))
c$$$      shp(4, 6)=xi*(-0.375+et*(-0.250+et*(0.125+ze*0.125)
c$$$     &  +ze*(-0.125+ze*0.125))+ze*(-0.250+ze*0.125)
c$$$     &  +xi*(0.125+et*(0.125+ze*0.125)+ze*0.125))
c$$$      shp(4, 7)=xi*(-0.375+et*(-0.250+et*(0.125-ze*0.125)
c$$$     &  +ze*(0.125+ze*0.125))+ze*(0.250+ze*0.125)
c$$$     &  +xi*(0.125+et*(0.125-0.125*ze)-0.125*ze))
c$$$      shp(4, 8)=-0.250+et*(et*(0.250-ze*0.250)
c$$$     &  +ze*(-0.250+ze*0.250))+ze*(ze*0.250)
c$$$     &  +xi*(-0.125+et*(-0.250+et*(-0.125+ze*0.125)
c$$$     &  +ze*(0.375-ze*0.125))+ze*(0.250-ze*0.125)
c$$$     &  +xi*(0.125+et*(0.125-0.125*ze)-0.125*ze))
c$$$      shp(4, 9)=xi*(0.500+et*(-0.500-0.500*ze)+ze*0.500
c$$$     &  +xi*(-0.250+et*(0.250+ze*0.250)-0.250*ze))
c$$$      shp(4,10)=xi*(0.250+et*(-0.250+ze*(ze*0.250))+ze*(-ze*0.250))
c$$$      shp(4,11)=xi*(0.500+et*(-0.500+ze*0.500)-0.500*ze
c$$$     &  +xi*(-0.250+et*(0.250-0.250*ze)+ze*0.250))
c$$$      shp(4,12)= 0.500+et*(-0.500+ze*(ze*0.500))+ze*(-ze*0.500)
c$$$     &  +xi*(-0.250+et*(0.250+ze*(-ze*0.250))+ze*(ze*0.250))
c$$$      shp(4,13)=xi*(0.500+et*(0.500+ze*0.500)+ze*0.500
c$$$     &  +xi*(-0.250+et*(-0.250-0.250*ze)-0.250*ze))
c$$$      shp(4,14)=xi*(0.250+et*(0.250+ze*(-ze*0.250))+ze*(-ze*0.250))
c$$$      shp(4,15)=xi*(0.500+et*(0.500-0.500*ze)-0.500*ze
c$$$     &  +xi*(-0.250+et*(-0.250+ze*0.250)+ze*0.250))
c$$$      shp(4,16)= 0.500+et*(0.500+ze*(-ze*0.500))+ze*(-ze*0.500)
c$$$     &  +xi*(-0.250+et*(-0.250+ze*(ze*0.250))+ze*(ze*0.250))
c$$$      shp(4,17)= 0.500+et*(et*(-0.500-ze*0.500))+ze*0.500
c$$$     &  +xi*(-0.250+et*(et*(0.250+ze*0.250))-0.250*ze)
c$$$      shp(4,18)=xi*(0.250+et*(et*(-0.250-ze*0.250))+ze*0.250)
c$$$      shp(4,19)=xi*(0.250+et*(et*(-0.250+ze*0.250))-0.250*ze)
c$$$      shp(4,20)= 0.500+et*(et*(-0.500+ze*0.500))-0.500*ze
c$$$     &  +xi*(-0.250+et*(et*(0.250-ze*0.250))+ze*0.250)
!
!     computation of the local derivative of the global coordinates
!     (xs)
!     This is the transpose of the Jacobian.
      do i=1,3
        do j=1,3
          xs(i,j)=0.d0
          do k=1,20
            xs(i,j)=xs(i,j)+xl(j,k)*shp(i,k) ! Note: this is transposed from shape20
          enddo
        enddo
      enddo
!
!     Normalize the rows of the transposed Jacobian
!
      do i=1,3
         row_norm = sqrt(xs(i,1)**2 + xs(i,2)**2 + xs(i,3)**2)
c         write(*,*) "CZM: row_norm = ", row_norm
         do j=1,3
            xs(i,j)=xs(i,j)/row_norm
         enddo
      enddo
!
!     computation of the normalized transposed jacobian determinant
!
      xsj=xs(1,1)*(xs(2,2)*xs(3,3)-xs(2,3)*xs(3,2))
     &   -xs(1,2)*(xs(2,1)*xs(3,3)-xs(2,3)*xs(3,1))
     &   +xs(1,3)*(xs(2,1)*xs(3,2)-xs(2,2)*xs(3,1))
!
!     (xsi) (inversion of xs)
!
      xsi(1,1)=(xs(2,2)*xs(3,3)-xs(3,2)*xs(2,3))/xsj
      xsi(1,2)=(xs(1,3)*xs(3,2)-xs(1,2)*xs(3,3))/xsj
      xsi(1,3)=(xs(1,2)*xs(2,3)-xs(2,2)*xs(1,3))/xsj
      xsi(2,1)=(xs(2,3)*xs(3,1)-xs(2,1)*xs(3,3))/xsj
      xsi(2,2)=(xs(1,1)*xs(3,3)-xs(3,1)*xs(1,3))/xsj
      xsi(2,3)=(xs(1,3)*xs(2,1)-xs(1,1)*xs(2,3))/xsj
      xsi(3,1)=(xs(2,1)*xs(3,2)-xs(3,1)*xs(2,2))/xsj
      xsi(3,2)=(xs(1,2)*xs(3,1)-xs(1,1)*xs(3,2))/xsj
      xsi(3,3)=(xs(1,1)*xs(2,2)-xs(2,1)*xs(1,2))/xsj
c$$$!
c$$$!     computation of the global derivatives of the shape functions
c$$$!
c$$$      do k=1,20
c$$$        do j=1,3
c$$$          sh(j)=shp(1,k)*xsi(1,j)+shp(2,k)*xsi(2,j)+shp(3,k)*xsi(3,j)
c$$$        enddo
c$$$        do j=1,3
c$$$          shp(j,k)=sh(j)
c$$$        enddo
c$$$      enddo
c$$$!
      xi = xi - 1               ! shift back to centered shape function

      return
      end
